Two-dimensional lattice-fluid model with water-like anomalies 
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We investigate a lattice-fluid model defined on a two-dimensional triangular lattice, with the 
aim of reproducing qualitatively some anomalous properties of water. Model molecules are of the 
"Mercedes Benz" type, i.e., they possess a Da (equilateral triangle) symmetry, with three bonding 
arms. Bond formation depends both on orientation and local density. We work out phase diagrams, 
response functions, and stability limits for the liquid phase, making use of a generalized first order 
approximation on a triangle cluster, whose accuracy is verified, in some cases, by Monte Carlo 
simulations. The phase diagram displays one ordered (solid) phase which is less dense than the 
liquid one. At fixed pressure the liquid phase response functions show the typical anomalous behavior 
observed in liquid water, while, in the supercooled region, a reentrant spinodal is observed. 

PACS numbers: 61.20.-p, 64.60.Cn, 64.60.My, 65.20.-l-w 
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I. INTRODUCTION 



Water is an anomalous fluid with respect to several 
thermodynamic properties [ll, |^ y| . At ordinary pres- 
sures the solid phase (ice) is less dense than the corre- 
sponding liquid, the liquid phase has a temperature of 
maximum density, while both isothermal compressibil- 
ity and isobaric heat capacity display a minimum as a 
function of temperature. Moreover, the heat capacity 
is unusually large. There is general agreement, among 
physicists, that an explanation of such anomalous prop- 
erties is to be found in the peculiar features of hydro- 
gen bonds, and the ability of water molecules to form 
such kind of bonds |3) I3- It is also widely believed 
that the same physics should be responsible of the un- 
usual properties of water as a solvent for apolar com- 
pounds [a, 13 1 tl^^.t is of the hydrophobic effect, of high 
importance in biophysics [g- Nevertheless, a compre- 
hensive theory which explains all of these phenomena 
has not been developed yet. A lot of work has been 
done in "realistic" simulations [3, llfl HB HJI i based on 
different interaction potentials, but they generally re- 
quire a large computational effort, and it is not always 
easy to understand which detail of the model is impor- 
tant to determine certain properties. On the contrary, 
simplified models generally need easier numerical calcu- 
lations and allow quite easily to trace connections be- 
tween microscopic interactions and macroscopic proper- 
ties El El in Hi mill El mm in in Hi. a sim- 
plified mechanism which has been proposed to describe 
the relevant physics of hydrogen bonding is the follow- 
ing one (see for instance Refs. 1^ 125,) . Hydrogen bond 
formation requires that the two involved molecules are 
in certain relative orientations and stay (on average) at 
a distance which is larger than the optimal distance for 
Van der Waals interaction. In other words there exists 
a competition between Van der Waals interaction (al- 
lowing higher density and higher orientational entropy, 
but resulting in a weaker bonding) and hydrogen bonding 
(requiring lower density and lower orientational entropy, 



but resulting in a stronger bonding). This simple mech- 
anism has been implemented in different models, both 
on- m m El 113 and off-lattice 13 , in 3 IH il as 
well as 2 dimensions |23, lla, |23. One of them is the 
2-dimensional Mercedes Benz model, originally proposed 
by Ben-Naim |14| , in which model molecules possess three 
bonding arms arranged as in the Mercedes Benz logo. In 
recent papers by Dill and coworkers |2J, |23, a similar 
(off-lattice) model has been simulated at constant pres- 
sure by a Monte Carlo method, allowing to describe in 
a qualitatively correct way several anomalous properties 
of liquid water and also of hydrophobic solvation. Never- 
theless, in view of investigations on the behavior of water 
in contact with other chemical species, as it happens for 
instance in several biological processes, it would be de- 
sirable to obtain an even simpler representation of the 
physics of hydrogen bonding. 

In this paper we investigate a model of the Mercedes 
Benz type on the triangular lattice, with a twofold pur- 
pose. As mentioned above, we are first meant to explore 
the possibility of obtaining a simpler model with the same 
underlying physical mechanism, and with qualitatively 
the same macroscopic properties. Moreover, we are inter- 
ested in extending the model analysis to the global phase 
diagram and in particular to the supercooled regime, in 
which water anomalies are thought to find an explana- 
tion. Such a detailed analysis is just made easier by in- 
creased simplicity. Working on a lattice, we have to resort 
to a trick to describe hydrogen bond weakening, when the 
two participating molecules are too close to each other. 
Such a trick is similar to the one proposed by Roberts and 
Debenedetti for their 3-dimensional model |2l |23 • The 
energy of any formed bond is increased (weakened bond) 
of some fraction by the presence of a third molecule on 
a site close to the bond (i.e., on the third site of the tri- 
angle). Due to the presence of only three bonding arms, 
it is not possible to distinguish between hydrogen bond 
donors and acceptors, but this seems to be of minor im- 
portance to the physics of hydrogen bonding |2J|. Let 
us notice that the model has the same bonding proper- 



ties as the early model proposed by Bell and Lavis |l3|. 
and the same weakening criterion as the model recently 
investigated by Patrykiejew and coworkers |2^ l27l|. but 
here non-bonding orientations are added. Such a feature 
is essential to describe directional selectivity of hydrogen 
bonds. 

The paper is organized as follows. In Sec. II we de- 
fine the model in detail and analyze its ground state. In 
Sec III we introduce the first order approximation in a 
cluster variational formulation, which we employ for the 
analysis. Sec. IV describes the results and Sec. V is de- 
voted to some concluding remarks. 



II. MODEL FORMULATION AND GROUND 
STATE 

The model is defined on a two dimensional triangular 
lattice. A lattice site can be empty or occupied by a 
molecule with three equivalent bonding arms separated 
by 27r/3 angles. Two nearest-neighbor molecules interact 
with an attractive energy — e (e > 0) representing Van 
der Waals forces. Moreover, if two arms are pointing to 
each other, an orientational term —rj (rj > 0) is added 
to mimic the formation of a hydrogen (H) bond. Due 
to the lattice symmetry, a particle can form three bonds 
at most and there are only 2 bonding orientations, when 
the arms are aligned with the lattice, while we assume 
that w non-bonding configurations exist {w is another 
input parameter of the model). Finally, the H bond en- 
ergy is weakened by a term cr]/2 (c e [0, 1]) when a third 
molecule is on a site near a formed bond. In the two 
dimensional triangular lattice there are two such weak- 
ening sites per bond, so that a fully weakened H bond 
energy turns out to be —(1 — cjrj. Let us notice that, in 
the above description, H bonding is a 3-body interaction. 
The hamiltonian of the system can be written as a sum 
over the triangles 
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where Tiijk is a contribution which will be referred to 
as triangle hamiltonian, and ir,iriiir" label site config- 
urations for the 3 vertices r,r' ,r" , respectively. Possible 
configurations are empty site (i = 0), site with a molecule 
in one of the 2 bonding orientations (i = 1, 2) or in one 
of the w non-bonding ones (i = 3) (see Tab. QJ. The 
triangle hamiltonian reads 

Hyfc = -einiTij + rijUk + rikrii) (2) 

-T][hij{l - crik) + hjk{l - cui) + hki{l ~ crij)], 

where rii is an occupation variable, defined as n^ = 
for i = (empty site) and n^ = 1 otherwise (occupied 
site), while hij = 1 if the pair configuration {i,j) forms 



a H bond, and hi 



otherwise. Let us notice that 



TABLE I: Possible site configurations, with corresponding 
labels (j) and multiplicities {wi). 
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say A, B, C, and i,j, k are assumed to denote configura- 
tions of sites placed on A, B, C sublattices respectively. 
Assuming also that A,B,C are ordered counterclockwise 
on up-pointing triangles (and then clockwise on down- 
pointing triangles), we can define hij = 1 if i = 1 and 
j — 2 and hij — otherwise. Let us notice that both 
Van der Waals {—enirij) and H bond energies {—rjhij), 
that are 2-body terms, are split between two triangles, 
whence the 1/2 prefactor in Eq. Q. On the contrary the 
3-body weakening terms {rihijcnk/2) are associated each 
one to a given triangle, and the 1/2 factor is absorbed 
in the prefactor. Let us denote the triangle configura- 
tion probability by Pijk, and assume that the probabil- 
ity distribution is equal for every triangle (no distinction 
between up- or down-pointing triangles) . Taking into ac- 
count that there are 2 triangles per site, we can write the 
following expression for the internal energy per site of an 
infinite lattice 
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The multiplicity for the triangle configuration (i,j,k) is 
given by WiWjWk, where Wi — w for i ~ 3 (non-bonding 
configuration) and Wi = 1 otherwise (bonding configura- 
tion or vacancy). 

Let us now have a look at the ground state properties of 
the model. In order to do so, let us investigate the zero 
temperature grand-canonical free energy lij° = u — ^p 
(/i being the chemical potential and p the density, i.e., 
the average site occupation probability), which can be 
formally written in the same way as the internal energy u 
of Eq. lO, by replacing the triangle hamiltonian Tiiju by 
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We find an infinitely dilute "gas" phase (G) with zero 
density and zero free energy, and an ordered "open 
ice" phase (Iq) with maximum number of H bonds per 
molecule. The latter configuration is realized through the 
formation of an open (honeycomb) H bond network with 
density 2/3 and free energy 



= — e — ?7 — 2/i/3. 



(5) 



triangle vertices are set on three triangular sublattices. 



Another possibility is the "closed ice" phase (Ic), in which 
all interstitial sites are occupied and all hydrogen bonds 



are fully weakened. The resulting free energy is 
ujI = -3e - 77(1 - c) - ^. 



(6) 



Let us notice that it is never possible to form 3 bonds 
in a triangle, which means that we have frustration. It 
is easy to show that the G phase is stable (wf > 0) for 
/i < /iG-i„, where 



MG-Io 



-3(6 + r,)/2, 



(7) 



the lo phase is stable (wf < and uj^ < u;^ ) for /^g-Io < 
/i < /ii^_i^, where 



Wo-ic = ^6e + 3c77, 



(8) 



and the Ic phase is stable (wj? < and wf < i^f ) for 
/i > /ii^_i^. The lo phase has actually a stability region, 
i.e., /XG-io < Mio-icJ provided 
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(9) 



which, in the worst case (c = 0), reads rj > 3e. We shall 
always work in the latter regime, which is the most sig- 
nificant one to describe real water properties. It is also 
possible to show that, at the transition point between the 
open and closed ice phases {fi = /xi^_i^), any configura- 
tion built up of a honeycomb H bond network with any 
number of occupied interstitial sites has the same free en- 
ergy. Hence we expect that the Iq — Ic transition does not 
exist at finite temperature, and actually we shall observe 
a unique ice (I) phase, in which the interstitial site oc- 
cupation probability gradually increases upon increasing 
the chemical potential. 

Let us finally notice that another possible phase is a ho- 
mogeneous and isotropic one in which the lattice is fully 
occupied and molecules can assume only bonding config- 
urations (i — 1,2). This "bonded liquid" phase, whose 
free energy coincides with that of the Ic phase in Eq. jnj, 
is observed in the u; = case, studied by Patrykiejew and 
others 26, 27]. In this scenario, non-bonding configura- 
tions are absent and the bonded liquid ground state has, 
for c 7^ 1, the same degeneracy as the Ising triangular 
antiferromagnet |23|. Nevertheless, in this work we shall 
deal with the case w ^ 1, which is relevant to describe 
H bond directionality. In this case the closed ice phase 
is entropically favored with respect to the bonded liquid 
phase, which cannot appear at finite temperature. In 
conclusion, because of the introduction of non-bonding 
configurations, the ground state degeneracy is removed 
at T = 0+, where only an infinitely dilute (gas) phase 
and a symmetry-broken (ice) phase are present. Such a 
phase behavior is closer to the one of water than the one 
obtained for w = 0. 



approximation on a triangle cluster, which we introduce 
in the framework of the cluster variation method. The 
cluster variation method is an improved mean-field the- 
ory based on an approximate expression for the entropy. 
In Kikuchi's original formulation 30] the entropy is ob- 
tained by an approximate counting of the number of mi- 
crostates. In a modern formulation 31] the approximate 
entropy can be viewed as a truncation of a cluster cu- 
mulant expansion. The truncation is justified by the ex- 
pected rapid vanishing of the cumulants upon increasing 
the cluster size, namely when the cluster size becomes 
larger than the correlation length of the system (the 
method necessarily fails near critical points) [S^]- The 
approximation is completely defined by the maximum 
clusters left in the truncated expansion, usually denoted 
as basic clusters. One obtains a free energy functional 
in the cluster probability distributions, to be minimized, 
according to the variational principle of statistical me- 
chanics. 

For our model we choose up-pointing triangles as ba- 
sic clusters (an analogous treatment works for down- 
pointing triangles). This approximation, which seems 
to be good in particular for frustrated models [SJ, ^M ' 
is easily shown to be equivalent to a first order approxi- 
mation on a triangle cluster ^^ . Let us notice that the 
internal energy is treated exactly, because the range of 
interactions does not exceed the basic cluster size, un- 
like the ordinary mean-field approximation. The grand- 
canonical free energy per site to = u — up — Ts (s being 
the entropy per site), can be written as a functional in 
the triangle probability distribution as 

333 
I3uj = '^'^'^w^WjWkPijk >^ (10) 

i=0 j=0 fc=0 

(SHi-jk + \npijk ~ 3 In {ptpfPk. 

where /3 = l/T (temperature is expressed in energy units, 
whence entropy in natural units) and pf- is the probabil- 
ity of the i configuration for a site on the X sublattice 
{X = A, B, C). The latter can be obtained as a marginal 
of the triangle configuration probability pijk, namely 

3 3 

j=0 k=Q 
3 3 

pf = ^^mwkPijk (11) 

•i=0 k=0 
3 3 

Pk = J2J2'^^^iP^3''- 
i=0 j=0 



III. FIRST ORDER APPROXIMATION 

We shall carry out the finite temperature analysis of 
the model mainly by means of a generalized first order 



The above expressions show that the only variational pa- 
rameter in Lu is the triangle probability distribution, that 
is the 64 variables {pijk}- 

The minimization of lu with respect to these variables. 



with the normalization constraint 

3 3 3 



^ ^ ^ W^WjWkPijk = 1, 
i=0 j=0 k=0 



(12) 



can be performed by the Lagrange niultipher method, 
yielding the equations 

p,,,=^^e-^^^-(p^fp^)^/^ (13) 

where ^, related to the Lagrange multiplier, is obtained 
by imposing the constraint Eq. p2l) : 



^ = J2J2J2 w.w.Wke-^^^^" {pfpfp'^) 



2/3 



(14) 
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Eq. (jl^ll is in a fixed point form, and can be solved numer- 
ically by simple iteration (natural iteration method [35J)- 
In our case the numerical procedure can be proved to 
lower the free energy at each iteration ;3.4 iS.Si] , and there- 
fore to converge to local minima. The solution of Eq. H13|) 
gives the equilibrium {pijk} values, from which one can 
compute the thermal average of every observable. Insert- 
ing these values into Eqs. (Q and (|10|l gives respectively 
the equilibrium internal energy and free energy. The lat- 
ter can be also easily expressed through the normaliza- 
tion constant as 



f3uj 



Ine, 



(15) 



whence ^ can be viewed as the approximate (single site) 
grand-canonical partition function. It is also worth men- 
tioning that Eq. (|13() preserves homogeneity (pf = pf ; 
Vi, X, Y), due to the invariance of Tiijk under cycle per- 
mutation of the subscripts (see Eqs. ||2J) and |0J). Let 
us finally notice that the free energy expression Eq. H10() 
can be also derived by considering the model on a tri- 
angular Husimi tree (triangle cactus) [3J| as a bulk free 
energy density, that is the free energy contribution far 
enough from the boundary, where an invariance condi- 
tion for the configuration probability of the triangles is 
assumed to hold. 



IV. RESULTS 

A. Phase diagrams 

In order to provide a first insight into the model, let 
us report in Fig. ^ the phase diagram in the chemical 
potential-temperature plane, for 7//e = 4, c = 0.5, and 
w = 50. Three phases can be observed: An ice (I) phase, 
with broken symmetry among the three sublattices, a 
liquid (L) phase and a gas (G) phase. The latter two 
phases preserve the sublattice symmetry but the liquid 
phase has a higher density. The ice phase has a lower 
density than the liquid phase, and its structure reminds 
that of ground state ice, with interstitial sites occupied by 
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FIG. 1: Temperature (T/e) vs. chemical potential (/^/e) 
phase diagram for w — 50, rj/e = 4, c = 0.5. G, L, and I 
denote the gas, liquid and solid (ice) phases respectively. CP 
denotes the critical point and TRP the triple point. 



molecules in non-bonding configurations. We can observe 
a triple point (TRP), in which the three phases coexist, 
and a gas-hquid critical point (CP). All displayed transi- 
tion lines are first-order. The above phase diagram shares 
several properties with the one of real water. Other crys- 
talline phases, such as a real close-packed ice, cannot be 
reproduced by the model. 

Let us now investigate the role of model parameters, 
by analyzing phase diagrams obtained for different val- 
ues. In Fig. I2K, v/^ f^nd c are left unchanged, while the 
number of non-bonding configurations w is varied within 
the interval [20, 100]. Upon increasing w, the liquid phase 
turns out to be more stable with respect to the ice phase, 
and the I-L transition temperature decreases. On the 
contrary, for lower w values, the I phase is increasingly 
stabilized and the I-L transition temperature increases. 
For w — 20 the whole L-G coexistence and also the crit- 
ical point disappears. Such a behavior can be explained 
by the fact that the L phase is characterized by a higher 
number of non-bonding molecules than the I phase, in 
which bonding molecules tend to form an ordered struc- 
ture. Therefore high w values largely increase the liquid 
phase entropy. 

In Fig. ISb, w and c are held fixed and the ratio rj/e 
is varied within the interval [3, 5]. Let us notice that we 
have restricted the investigation to cases in which the 
orientational (H bond) interaction is stronger than the 
non-orientational one, which is the case for real water. 
It turns out that the ratio rj/e affects the stability of the 
I phase with respect to both the G and L phases. In 
fact higher values of rj means stronger H bond, which fa- 
vors the I phase, that is the only extensively H-bonded 
phase. On the contrary the L and G phases are dom- 
inated by non-oriented interactions with coupling con- 
stants e, therefore both these two phases are unfavored 
by high rj/e values. Even in this case the L-G coexistence 
may become metastable. 

The ice phase at high pressures has maximum density 



and number of weakening molecules per H bond. Raising 
c, the stability of this configuration is lowered with re- 
spect to the liquid phase with few H bonds. This is shown 
in Fig. 13; where rj/e and w are fixed and the weakening 
parameter c is varied in its interval of definition [0,1]. 
This trend is reversed for low w values (w = as well), 
because in the latter case the liquid has the maximum 
number of fully weakened bonds. 

In the next part of this work we focus on a par- 
ticular choice of parameters (w — 20, rj/e — 3 and 
c = 0.8) which, from the above analysis, turn out to 
correspond to a water-like phase diagram. Fig. |31 shows 
the temperature-pressure phase diagram, and Fig. 0| the 
temperature-density phase diagram. Let us notice that 
pressure P is simply given by P = —uj (the volume per 
site is assumed to be equal to 1, i.e., pressure is expressed 
in energy units), due to the fact that the free energy has 
been defined as a grand-canonical potential. 
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FIG. 3: Pressure (-P/e) vs. temperature (T/e) phase diagram 
for w = 20, TTj/e — 3 and c — 0.8. Solid lines denote first 
order transitions, a dashed line denotes the TMD locus, and 
a dash-dotted line denotes the stability limit for the liquid 
phase. The inset displays, in addition, the locus of divergence 
of the density response functions at low temperature (solid 
line) with its "critical" point and the Kauzmann line (dashed 
line). 
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FIG. 2: The same phase diagram as in Fig. (dashed lines) 
compared to different parameter choices: (a) lo = 20 (solid 
lines) and w = 100 (dash-dotted lines); (b) rj/e — 5 (solid 
lines) and rj/e = 3 (dash-dotted lines); (c) c — 0.8 (solid 
lines) and c — 0.2 (dash-dotted Unes). 
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FIG. 4: Temperature (T/e) vs. density (p) phase diagram 
for w = 20, rj/e = 3 and c = 0.8. Solid lines denote phase 
boundaries; a thin dashed line corresponds to the triple point. 
Phase labels as in Fig. double labels denote two-phase co- 
existence regions. 



B. TMD locus and stability limits 

One of the water anomalies that the present model is 
able to reproduce is the temperature of maximum den- 
sity (TMD) along isobars for the liquid phase. Joining 
TMD at different pressures defines the so called TMD 
locus, which is a negatively sloped line in the T-P phase 
diagram of real water. We determine the TMD locus nu- 
merically, by adjusting the chemical potential in order 
to fix the pressure and then imposing that the (isobaric) 
thermal expansion coefficient vanishes. 

The limit of stability of the liquid phase (spinodal) is 
the locus in which the metastable liquid ceases to be a 



minimum of the free energy, and becomes a saddle point. 
The stabihty hmit can be obtained by studying the eigen- 
values of the hessian matrix of the free energy |3fil | 
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Let us notice that, when the liquid phase stability is lost 
(some eigenvalue of the above matrix vanishes), also the 
corresponding fixed point of the natural iteration equa- 
tions (|13|l becomes unstable. In order to determine the 
stability limit with respect to the symmetry-broken ice 
phase, it is sufficient to impose homogeneity during the 
iterative procedure, which is done by replacing Eqs. 111|) 
with 



pf=P? 



3 3 

j=0 k=0 



Pi]k + Pki] + Pjki 



(17) 



This trick cannot be applied when the liquid stability is 
lost with respect to a homogeneous phase, because the 
liquid fixed point of equations (|13|l becomes definitely 
unstable, due to divergence of the density response func- 
tions. In the latter case the spinodal is determined by 
solving the eigenvalue problem for the hessian matrix 
rewritten by forcing the homogeneity condition (|17|l . 

The results are shown in Fig. 13 The stability limit 
of the liquid with respect to the gas phase starts from 
the critical point and reaches a minimum in the negative 
pressure region. After this point the line becomes neg- 
atively sloped and joins continuously the stability limit 
with respect to the ordered ice phase. The TMD locus 
intersects the limit of stability in its minimum in the 
T-P plane, according to the predictions of Speedy and 
Debenedetti il @ iHl ES IS El El , based on thermo- 
dynamic consistency arguments. In fact the TMD locus 
causes the liquid limit of stability line to retrace, giving 
rise to a tensile strength maximum and to a continu- 
ous boundary. Let us recall that, while at the stability 
limit with respect to the gas phase, the density response 
functions diverge, this is not the case at the stability 
limit with respect to the ordered phase. Nevertheless we 
can observe that the density response functions tend to 
diverge also upon decreasing temperature, as observed 
experimentally. The locus of divergence, terminating at 
some kind of critical point, can be defined, in the frame- 
work of a simplified variational free energy forced to de- 
scribe a homogeneous system, as an additional stability 
limit with respect to a low density liquid phase. Such 
"phase" corresponds to a saddle point of the original 
(not symmetrized) free energy, unstable with respect to 
the solid phase. As the low pressure solid phase reminds 
the ground state "open ice" structure, which is three-fold 
degenerate, the triangle probability distribution of the 
low density liquid phase turns out to be essentially an 
arithmetic average over the three ice distributions. The 
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FIG. 5: Response functions at constant pressure (P/e = 1) 
as a function of temperature {T/e): (a) thermal expansion 
coefScient {eap); (b) specific heat {cp); (c) isothermal com- 
pressibility (ekt). 



unphysical nature of this solution is also reflected in its 
negative entropy. The divergence locus, together with the 
locus at which the liquid phase entropy vanishes (Kauz- 
mann line), are shown for completeness in the inset of 
Fig. 13 Upon increasing temperature the divergence lo- 
cus meets the spinodal tangentially and they become the 
same curve ending in the "true" gas-liquid critical point. 



C. Response functions 

Let us now investigate the density response functions 
and the specific heat of the liquid at constant pressure 
P/e = 1 (pressure is kept fixed by numerically adjusting 
the chemical potential /x). It turns out that these func- 
tions display an anomalous behavior similar to that of 
real liquid water. The first response function we consider 
is the thermal expansion coefficient ap = {—d\np/dT)p, 
which is proportional to the entropy-specific volume 
cross-correlation. For a typical fiuid ap is always positive 
because if in a region of the system the specific volume 
is a little larger then the average, then the local entropy 
is also larger, i.e., the two quantities are positively cor- 
related. On the contrary, for our model ap (Fig. [SJl) 
displays an anomalous behavior. As temperature is low- 
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FIG. 6; Gas-liquid transition at fixed temperature {T/e — 
1.05), upon varying the chemical potential (l-i/e): First order 
approximation results (solid line) compared to Monte Carlo 
simulations (scatters), for w — 20, -q/e = 3, and c — 0.8. The 
inset displays the Binder cumulant minimum, together with 
the transition point predicted by the first order approximation 
(vertical line). 
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FIG. 7: Stability limits from first order approximation (thick 
lines) and homogeneous nucleation points from Monte Carlo 
simulations (scatters), for w — 20, rj/e — 3, and c — 0.8. 
Open circles and a solid line denote the stability limit to the 
ice phase, filled circles and a dashed line the stability limit 
to the gas phase. Thin dash-dotted lines denote equilibrium 
phase boundaries. 



ered ap vanishes (at the TMD), becomes negative, and 
finally tends to diverge. As previously mentioned, diver- 
gence can be observed only for pressure values less than 
some "critical" pressure. Anyway, before divergence is 
actually reached, the liquid loses stability with respect 
to the ice phase. 

The trend of the isothermal compressibility kt = 
{d\np/dP)T is also anomalous (Fig. |S|:). For a typical 
liquid Kt decreases as one lowers temperature, because 
it is proportional to density fluctuations, which decrease 
upon decreasing temperature. On the contrary, in Fig.lsj; 
we can observe that kt, once reached a minimum, be- 
gins to increase upon decreasing temperature. Such a 
behavior is observed in real liquid water. An analogous 
behavior characterizes the constant pressure specific heat 
cp = {-Td^ix/dT^)p (Fig.Eb). 



D. Numerical simulation 

We have studied the model in the first order approxi- 
mation to obtain easily detailed information about phase 
diagrams and in particular the metastable region. In 
order to check this approximation and obtain an esti- 
mate of its quantitative accuracy, we have also performed 
some (grand-canonical) Monte Carlo simulations on a 
60 X 60 triangular lattice with periodic boundary con- 
ditions. From the very beginning, we have chosen quite 
a low number of non-bonding configurations for our anal- 
ysis {w = 20), in order to increase the speed of simula- 
tion dynamics. In fact a lower w value corresponds to a 
smaller configuration space. We report some results in 
the following. 

In Fig. we show a first order transition between the 
gas and the liquid phases along a constant temperature 



path T/e = 1.05, quite less than the critical temperature. 
At the critical point the correlation length increases and 
the approximation may give worse predictions. Fig. El 
suggests that the first order approximation well local- 
izes the transition and that far enough from the critical 
point its predictions are nearly quantitative. Of course 
Monte Carlo simulations display smooth density varia- 
tions, due to finite size effects, but the Binder cumulant 
(inset), displaying a minimum, gives evidence of a first 
order transition. 

The reentrance of the liquid stability limit, which is one 
of the striking features of the (metastable) phase diagram 
of this model, is also confirmed by simulations. Perform- 
ing simulations in the metastable region, the spinodal 
has been determined by an arbitrary criterion for the life 
time of the metastable phase (100 Monte Carlo steps), as 
it has been done in previous studies 22] . Such a criterion 
allows us to find the kinetically controlled limit of super- 
cooling (homogeneous nucleation locus), shown in Fig.|3 
along with the corresponding first order approximation 
result. Both methods show a reentrant spinodal form- 
ing a continuous boundary. The simulations also confirm 
the distinction between liquid limit of stability with re- 
spect to the gas or to the ice phase, as in the first order 
approximation. 



V. DISCUSSION AND CONCLUSIONS 

In this paper we have investigated a 2-dimcnsional lat- 
tice model in which model molecules possess three equiv- 
alent bonding arms, and bonding energy depends on the 
presence of neighbor molecules, giving rise to a 3-particle 
interaction. The observed behavior is qualitatively sim- 
ilar to that of water, exhibiting the correct anomalies. 
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Upon supercooling, kt and cp increase and ap becomes 
negative and large in magnitude. Nevertheless, at ordi- 
nary pressures (less than the critical pressure) the density 
anomaly {ap — 0) is found in the metastable liquid re- 
gion. We have also determined the spinodal limits to the 
liquid state, and pointed out the relationship between 
these limits and the TMD locus. The growth in the re- 
sponse functions upon decreasing temperature can be in- 
terpreted on the basis of a reentrant spinodal scenario. 
The liquid-gas spinodal meets the TMD locus at the reen- 
trance point, as required by thermodynamic consistency. 
Actually the reentrant spinodal conjecture is one of the 
possible theoretical explanations of water anomalies, and 
some experimental results are consistent with this expla- 
nation |4J|. Nevertheless, it is important to note that, for 
the specific case of water, alternative interpretations of 
the stability problem exist, based on the second critical 
point conjecture 4]. The latter, supported by molecular 
dynamics simulations 12], seems to be more consistent 
with the existence, in the negative pressure region, of a 
monotonic liquid-gas spinodal and a reentrant TMD lo- 
cus. On the contrary, our model displays a metastable 
liquid state which is bounded by a spinodal both at pos- 
itive as well as negative pressures, forming a continuous 
boundary. The lower temperature part of the bound- 
ary is the limit of stability with respect to the ordered 
ice phase, while the higher temperature part is the limit 
of stability with respect to the gas phase. While the 
response functions diverge at the liquid-gas spinodal, at 
the liquid-solid spinodal they do not, even if they tend to 
higher values. Anyway, in our framework, it is also pos- 
sible to investigate the behavior of the unstable liquid 
(a saddle point of the variational free energy) and de- 
termine the locus of divergence. The latter always turns 
out to lie at a temperature less than the limit of stabil- 
ity, according to experiments [43 ■ It also turns out that 
the divergence locus terminates at some kind of critical 
point, meaning that response functions should not show 
divergent-like behavior for pressure values greater than 
some critical pressure. 

Let us notice that a previous lattice model on the 3- 



dimensional body centered cubic lattice had pointed out 
a qualitatively similar behavior ^^- Nevertheless, in 
such a model, orientational degrees of freedom of water 
are not treated explicitly and two equivalent sublattices 
are artificially distinguished by the hamiltonian. This is 
necessary to favor an open structured phase. Moreover, 
the analytical treatment is based on the determination of 
a temperature dependent 2-particle interaction. On the 
contrary in our model there exists an explicit, though 
simplified, modelling of hydrogen bonding and no tem- 
perature dependent interaction is introduced. The open 
structured phase is favored in principle by the triangular 
lattice structure. 

We have mentioned in the Introduction that the 
present model is actually an extension over an early 
model proposed by Bell and Lavis |l3l | (corresponding 
to the case in which w = and c = 0) and over a 
recent model investigated by Patrykiejew and cowork- 
ers |2^|23 (corresponding to w = 0). The former model 
in the same approximation actually displays, for 77/e > 3, 
a density anomaly (without singularities), but we have 
verified that the anomaly occurs in a negative entropy 
region. The latter model shows an unrealistic phase di- 
agram, in which, for high enough pressure, the liquid 
phase extends its stability region down to zero tempera- 
ture. In the present work we have shown that the addi- 
tion of non-bonding configurations to such a simple class 
of 2-dimensional lattice models allows us to reproduce a 
qualitatively correct water-like behavior. Moreover, this 
result has been obtained in a computationally much sim- 
pler way than a conceptually similar model with con- 
tinuous degrees of freedom, that is the Mercedes-Benz 
one. The latter model is highly appealing, because of its 
ability to explain most phenomena related to hydropho- 
bicity ^3- Therefore it would be interesting to analyze 
also the properties of the present model for a solution 
of an inert (apolar) solute, whose peculiar properties are 
thought to be strictly related to hydrogen bonding. This 
goes beyond the scope of the present paper and will be 
the subject of a forthcoming article. 
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